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This paper studies various Hopf bifurcations in the two-dimensional plane Poiseuille problem. 
For several values of the wavenumber a, we obtain the branch of periodic flows which are born at 
the Hopf bifurcation of the laminar flow. It is known that, taking a « 1, the branch of periodic 
solutions has several Hopf bifurcations to quasi-periodic orbits. For the first bifurcation, previous 
calculations seem to indicate that the bifurcating quasi-periodic flows are stable and go backwards 
with respect to the Reynolds number, Re. By improving the precision of previous works we find that 
the bifurcating flows are unstable and go forward with respect to Re. We have also analysed the 
second Hopf bifurcation of periodic orbits for several a, to find again quasi-periodic solutions with 
increasing Re. In this case the bifurcated solutions are stable to superharmonic disturbances for Re 
up to another new Hopf bifurcation to a family of stable 3-tori. The proposed numerical scheme 
is based on a full numerical integration of the Navier-Stokes equations, together with a division by 
3 of their total dimension, and the use of a pseudo-Newton method on suitable Poincare sections. 
The most intensive part of the computations has been performed in parallel. We believe that this 
methodology can also be applied to similar problems. 



I. INTRODUCTION 



The theory of hydrodynamic stability is one of the main topics in fluid mechanics. Poiseuille as well as Taylor- 
Couette flow are test problems where it is possible the evaluation of different analytical and numerical methods, 
due essentially to the simplicity of their geometry. The dynamics of plane Poiseuille flow departs from the laminar 
flow. The stability of the laminar solution to infinitesimal disturbances has been analysed linearly and gives rise to 
the Orr-Sommerfeld equation. This equation has been studied by several authors as Thomas [T], Orszag [5], and 
Maslowe [3] among others, and it is well understood. The critical Reynolds number of the linear theory, Rccr = 
5772.22 for the wavenumber a = 1.02056, has been obtained by this approach. However, as experiments of |Carlson7 



Widnall, and Peeters| |4j , Nishioka and Asai [5] , and Alavyoon, Henningson, and Alfredssonj [6] showed, transition 



to turbulence is observed for Reynolds number 1000, what motivates that finite- amplitude disturbances originate 
the transition. The understanding of the transition to turbulence has been conjectured by SafFman [7] to depend on 
intermediate vortical states and turbulence takes place due to their three-dimensional instability. In recent years, 
authors have also payed attention to subcritical transition models based on transient optimal growth (see Schmid and 
Henningson [5], for instance). Examples of vortical states are periodic[l5] flows in time or space, among which can be 
mentioned: two-dimensional travelling waves, secondary flows in two or three dimensions (for them the flow rate and 
the pressure gradient are constants) and quasi-periodic solutions. Ehrenstein and Koch [9 discovered a new family 
of secondary bifurcation branches in dimension 3, which contains only even spanwise Fourier modes and reduces the 
critical Reynolds number (defined in terms of the averaged velocity across the channel) to Req^^ ~ 1000 as observed 
in experiments. 

Two-dimensional disordered motion is associated with the large scales of some turbulent flows, so there probably 
exist attractors for those two-dimensional flows. Besides, two- and three-dimensional states can compete and coexist 
in the final flow (cf. Jimenez [TD] and the references therein). In spite of the fact that transition to turbulence is a 
three-dimensional phenomenon, there are many properties of the two-dimensional flows observed in fully turbulent 
three-dimensional flows such as wall sweeps, ejections, intermittency and bursting, as Jimenez [IF showed. The two- 
dimensional case has attracted the attention of many authors but it is not completely understood as the problem of 
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two-dimensional transition to turbulence proves. Due to Squire s |12| theorem, to every three-dimensional perturbation 



of the linearized Navier-Stokes equations for a given Re, a, it corresponds a two-dimensional one for some a ^ a 
and Re ^ Re, so the critical Re for the linear theory must be attained by a two-dimensional flow. This result has 
been one of the main reasons to firstly try to understand the two-dimensional case, apart from the obvious easiness 
of computations compared to the three-dimensional situation. In addition, some of the properties obtained from the 
two-dimensional case can also provide new insight for three-dimensional flows. 

In this work we intend to analyse the dynamics of an easily treatable problem without domain complexities as 
is the case of the two-dimensional plane Poiseuille flow. Different levels of bifurcation to respective vortical states 
are considered, starting at the basic parabolic flow. From it, a family of travelling waves is born subcritically (see 



IV C I for a ~ 1. There are many papers concerning this kind of waves: Soibelman and Meiron fl3' gave an excellent 
review about it. As a starting point for our computations we have also reproduced the calculations to find the 
travelling waves for several values of a. Jimenez jTUl [H] and Soibelman and Meiron [13] obtained the next level 
of bifurcation to quasi-periodic solutions. Employing full numerical simulation in time, Jimenez |10L lllj computed 
different attractor flows with a moderate number of Chebyshev and Fourier modes. On the other hand, Soibelman 
and Meiron [13 implemented an algebraic approach to capture stable and unstable quasi-periodic flows, but the 
number of modes used were not enough to give good results and they were not able to carry out the stability analysis. 
The method implemented in the present work combines both: we solve a stationary problem to compute travelling 
waves for an observer moving at an appropriate speed, whereas the quasi-periodic flows are found by means of full 
numerical integration of the Navier-Stokes equations. Through algebraic manipulations, we express the discretized 
Navier-Stokes system only in terms of the stream component of the velocity. As a consequence, the dimension of 
the system is divided by 3, reducing considerably the computational effort. Using the numerical integrator, we have 
built a Poincare section of the flow, in order to apply a pseudo-Newton method for obtaining also unstable quasi- 
periodic solutions. These unstable intermediate states of the flow provide a highly useful insight into the transition 
process, as exemplified by secondary bifurcations in shear flows (see Casas and Jorba for instance). The spatio- 
temporal symmetries of the channel allows the reduction of quasi-periodic flows with two-frequencies to periodic flows 
in the appropriate Galilean reference. The quasi-periodic solutions found in this work correspond to the first two 
Hopf bifurcations of travelling waves for the case of constant pressure drop through the channel, and the first Hopf 
bifurcation when the mass flux is held constant. The property of behaving as time-periodic flows if we take a suitable 
Galilean reference, simplifies enormously the search of this kind of solutions. For them, the associated return time 
to the Poincare section is roughly 10000 time units at the first Hopf bifurcation for constant pressure, what makes 
the temporal integration very costly. The considered numerical procedure utilizes a parallel algorithm to evaluate the 
different columns of a Jacobian matrix, needed in the application of pseudo-Newton's method for the continuation 
of quasi-periodic solutions. We find that on the analysed Hopf bifurcations for both constant pressure and constant 
flux formulations, there exist quasi-periodic flows with increasing Re for some range of a and with decreasing Re 
for some other a: the bifurcations are supercritical or subcritical respectively. On the first bifurcation for constant 
pressure, we have traversed a curve of unstable quasi-periodic solutions. On the remaining bifurcations, there are 
stable quasi-periodic solutions to disturbances with the same wavenumber a and likewise, for Re sufficiently large, we 
have obtained unstable solutions. 

Once we have situated the different studies concerning Poiseuille flow, in the next section we pose the concrete 
terms that define the plane Poiseuille problem in two dimensions, together with their equations for both cases of 
constant pressure and flux. Next in§ |III| we explain the main details of the numerical methods. In § |IV| we review some 
results of the Orr-Sommerfeld equation and obtain, for several values of a, the bifurcating solutions of time-periodic 
flows. From these we analyse in§[V]the bifurcating branches to quasi-periodic solutions at the above-mentioned Hopf 
bifurcations. Finally in § |VI| we point out some conclusions. 

II. POISEUILLE FLOW 

We consider the flow of a viscous incompressible two-dimensional fluid, in a channel between two parallel walls, 
governed by the Navier-Stokes equations together with the incompressibility condition 

^ + (u.V)u = -Vp-f — Au, V-u = 0, (I) 

where u = u{x,y,t) = {u,v){x,y,t) represents the two-dimensional velocity, p = p{x,y,t) the pressure and Re the 
Reynolds number. As boundary conditions we suppose no-slip on the channel walls at y — ±1 and, at artificial 
boundaries in the stream direction period L, i.e. 



uix,±l,t) = v{x,±l,t) = 
(u,v,pl)[x -f L,y,t) = (u,v,p'){x,y,t) 



xeR, ye [-1,1], t^O, (2) 



3 



being p' = p + Gx, for G = G{t) the mean pressure gradient on the channel length, L, in the streamwise direction. 
For the system previously described there is a time-independent solution known as the basic or laminar flow that has 
a parabolic profile, namely 



ubiy) = 1 - 



0, Vpb 



Magnitudes in ([T])-([2]) are non-dimensional. We consider the two typical formulations used to drive the fluid: 
fixing the total flux Q, or the mean pressure gradient G, through the channel. For each of them we obtain a different 
definition of Re ~ hUc/v namely, Req — "iQ/Av and Rcp = Gh^ /2pv^ respectively, where, in dimensional magnitudes, 
h represents half of the channel height, Uc the velocity of the laminar flow in the centre of the channel, and v and p 
the constant kinematic viscosity and density. For a given laminar flow, i.e. letting Uc fixed, both definitions of the 
Reynolds number coincides with Re = hUc/y- That is not the case for secondary flows, defined as the ones having 
constant flux and mean pressure gradient through the channel. If we consider such a flow u(a;, y), expressed for each 
formulation by means of respective Fourier series 

uQ(x,y) =^u^(y)e''="% uP(x,2/) = ^ u^y)e"="^ 

fcez fcgz 

then, using the notation [/]^ :— f{h) — ,f{a), it is easy to check that (see for instance Casas (TS]) 

1 



Rep 
Req 



dvl 
dy 



J -1 



Reg 



and the corresponding relationships between velocities and pressures 



uPix,y) 



Re, 



u'3(a;,y), pP{x,y) 



vP{x,y) dy, 



(3) 



(4) 



We will employ later that periodic conditions at artificial boundaries in the stream direction, yield a great simpli- 
fication in the structure of the flow: quasi-periodic solutions may be viewed as periodic flows, and periodic solutions 
as stationary ones, if the observer moves at adequate speed c, in the stream direction. For this reason we perform the 
change of variable x = x — ct, which (writing again x instead of x) turns system ([T]) into: 



du 
dt 

dv 
~di 



du 


du 
dy 


dp 


dx 


dx 


dv 


dv 
dy 


dp 


dx 


dy 


du 
dx 


dv 
dy 


0, 



d^u d^u 
dx"^ dy"^ 

d'^v d'^v 



(5) 



together with boundary conditions as in (H). We can recover ([T]) by simply taking c = in ([s]) . 



III. NUMERICAL APPROACH 



Let us now describe the numerical procedure. For system ([5| we want to follow the temporal evolution of an 
initial flow subjected to the incompressibility condition, V • u = 0, and boundary conditions ([2|. To this end we 
use a spectral method to approximate velocities u,v and pressure deviation p', which from now on we consider 
non-dimensional quantities. We recall that p = p' — Gx and as it is easily obtained (see for example Casas P3^) 



G = -- 



1 

2R^c 



duo 
. dy _ 



G = 



2 

Rer. 



(6) 



respectively for the constant flux or pressure cases, so in the first one the mean pressure gradient varies with time 
and it is constant for the second one. 
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Spatial discretization. We use a standard Fourier- Galerkin, Chebyshev-collocation approach (cf. Canute et al. 
[16j ) in order to discretize x^y derivatives. In this way, we consider Fourier series (with a = 2t: / L the parameter 
wavenumber) : 



N 



iu,v,p'){x,y,t) = ^ {uk,Vk,Pk)iy,t)e'''°"', xeR, ye[-l,l], i > 0, 

k=-N 



which substituted in ([5| gives rise to a system of partial differential equations for the Fourier coefficients {uk,Vk,Pk): 

d^Uk ' 



duk 


' , , du 


du 
dy 
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^ Re 


dvk 
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dv 
dy 




_ dpk ^ 
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Re \ 




ifcauj. 


^ dvk 
dy 


= 0, 





, 2 2 » , 9^i'k 
-k a Vk + -TT-iT- 
dy^ 



(7) 



where —N ^ k ^ N, [-J^. stands for the order fcth Fourier coefficient of [•], Sqq = 1, and Sko = for fc 7^ 0. Because 
u,v,p' are supposed to be real functions, it is enough to consider modes Uk,Vk,Pk for k — 0,...,N in ([T]). The 
corresponding no slip boundary conditions in ^ are now written as 



iuk,Vk){±l,t) ^0, for i ^ and A: = 0, 



(8) 



The previous system is imposed at two different sets of Chebyshev abscissas to avoid indeterminacy, namely ym — 
cos(7rm/M) (velocities and momentum) for to = I,...,M — I, and ym+1/2 — cos(7r(m + l/2)/M) (pressure and 
continuity) for m = 0, . . . , M — I. 

Reduced equations. To emphasize the linear character of some operations, we now write system ([t]) as 



Uk = 



Vk 



du du 

^"'-'^dx+'^ay 

. dv dv 



D,,kC^'C2Pk + ^(^xfe + C^'DlOuk + 4oG, 



C^'DyC2Pk + -^^{Dlk + C^'DlC^)vu, 



D^kC^^CiUk + C^^DyCiVk = 0, 



(9a) 

(9b) 
(9c) 



for k — 0, . . . , N , where we have taken out of [■]f.,Uk,Vk,Pk for convenience. In ([9|) we have represented vectors of 
values Uk,Vk at the grid ym and pk at the grid ym+1/2', Ci,C2 are the corresponding matrices of cosines transforms 
for grids ym and ym+1/2, and D^k, Dy denote the respective matrices of partial derivatives in x, y. 

From (9c| we obtain a matrix Tk that carries out the transformation Vk = T^Uk where Uk = (uk,i, ■ ■ ■ ,Uk^M-2Y 
and Vk = (Uk,M-i,Vk^i, ■ ■ ■ ,Vk,M-iY ior k — 1, . . . , N . For fc = 0, from the continuity equation in ([t]), we obtain 
dva/dy — 0. Applying boundary conditions, vq{±1) — 0, we get VQ{y) = 0. This implies Wo,i = ■ ■ ■ — Wo,7\/-i = 0- 



For k — I, ... ,N we introduce the notation 

Uk = 

Vk = 



. .du du 



dv dv 

[u — c)- — h V 



1 

Re 
1 



{Dlk 



C^^DlCi)uk + SkoG, 



+ ^iDik+C^'D'yC,)vk, 



dx dy ^, Re 

Uk = {Uk){l^...^M-2}, Qk — {DxkCi ^C2){l,...3,/-2}, 



Vk = 



{Uk){M-l} 

Vk 



n _ ( {DxkC-^ ^C2){M-i} 



where A^^^ ^^j stands for rows ii, . . . ,i„ of matrix A. Equations (9al and (9b) can be now expressed as 



Uk = Uk - QkPk, 
Vk = Vk- QkPk- 
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The matrix Qk turns out to be an M x M invertible matrix. Consequently, from the second equation we obtain 
Pk = Qk^{Vk — Vk), which substituted into the first one yields 

uk = tJk - QkQk\Vk - ik) - Uk - QkQk\Vk ~ TkUk). 

Finally letting Pk = QkQZ^ , we can also invert / — PkTk, and thus we may solve for Uk 



Uq = Uq, 

nk^{I-PkTk)'HUk-PkVk), fc = l, 



(10) 



where / is the identity matrix of dimension M 

mind the substitution Vk = TkUk, we observe that system (10) does not depend on Vk nor pk 
and Uk for k 



2 and we have extended the definition of Uk for k ~ 0. Bearing in 

it only depends on uq 



1, . . . , N. In addition, due to the elimination of pressure in ( 10 1, we avoid the indeterminacy caused by 



an additive constant. However this indeterminacy has no effect upon the pressure gradient. Likewise this formulation 
saves the problems in the imposition of consistent initial conditions with the incompressibility. At the same time the 
stability analysis is simplified from (10 1. 



Temporal evolution. Once removed v and p from ([9|, in ( 10 ) it just remains to discretize temporal derivatives. We 
can express ( 10 1 as 



Uk 



^ Ck[uk) +Mk{uQ,. . . ,un), fc = 0, 



(11) 



where uo = and Ck, A4 corresponds respectively to linear and nonlinear terms in uqj ■ • • t^n on the right hand 
side of (10 1. We adopt a usual scheme for advection-diffusion problems: letting -uJJ be Uk at the time instant n/S.t for 



some fix time step Ai, we approximate Nl 
by an implicit one (Crank-Nicolson), so that (11) yields 



-■Mk{i 



''N 



) by an explicit method (Adams-Bashforth) and Ck{u]t) 



^k - ^'^fcl^fc ) - 



At 



(12) 



For the kind of solutions treated in this work and moderate values of Re < 10000, we have verified local errors 
originated in ( 12 ) from the time discretization. For that purpose, we approximate temporal derivatives by central 
finite differences and then improve precision by means of extrapolations. In all tested cases we have found errors 



0((Ai)^), which is in agreement with the discretization errors in (12). For some flows considered in §|v]it has been 
necessary to reduce Ai to avoid overflows in u(t). 

We apply system (12) to the two formulations described in §|llj namely, constant flux and constant mean pressure 
gradient. The imposition of constant flux Q = 4/3 (a linear condition) allows us to reduce by one the number of 
unknowns in ug. Therefore the number of equations is also reduced by one. This condition is related to the formula 
derived for G in ([6|, which depends linearly on uq and thus it is included in Lq. On the other hand, in the constant 
pressure case, the value of G is held constant and so it is a nonlinear term. Taking into account that in (10), uo 
has only real components but for fc = 1, . . . ,iV, it is a complex vector, we conclude that the block for A: = has 
dimension M — 2 or M — 1 respectively for Req and i?ep formulations, and dimension M — 2 for the 27V remaining 
real blocks. In summary, each time step, ( 12 ) implies the solution of a block diagonal linear system of total real 
dimension (2A^ + 1)(A/ — 2) in the constant flux case and (2iV + l)(Af — 2) + 1 in the constant pressure one. That 
means a rough division by 3 in the dimension of the whole system ([7| . In what follows we denote a solution flow at 
time t as U{t) = {uq, ... ,UN){t) G for K = {2N + 1)(M - 2) + 1 or = (27V + l)(Af - 2), according to the two 
above-mentioned cases. 



IV. PERIODIC SOLUTIONS 
A. The Orr— Sommerfeld equation 



Before applying the previously described numerical scheme, we make some considerations about the linearized 
stability of the laminar flow and time-periodic solutions. We start from the linearization of the vorticity equation 
around the basic flow, which is known as the Orr-Sommerfeld equation 

[ub - -)(</)" - a2</,) - u'^^p = ^(0(4) _ 2a^(b" + a^)- (13) 
a laRe 
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FIG. 1: Neutral stability curve (in blue) for the laminar solution, using n < 1000 discretization points. For each pair (Re, a) 
in this curve, the most unstable eigenvalue A is purely imaginary. The curve splits the Re-a plane in two stability regions as 
shown in the graph: the green one is stable and the red one unstable. 



It is a fourth order ordinary differential equation on (j) — 4>{y) as eigenfunction, with A as eigenvalue, and boundary 
conditions 0(±1) = (/''(±1) = 0. For each Re and a, (131 represents an eigenvalue problem on A and (j). In this way 
if A = Ar + iAi is a complex eigenvalue with A^ > 0, then the laminar flow is unstable to infinitesimal disturbances 
according to the linear theory. 

We have employed finite differences to approximate 4>{y) and its derivatives in an uniform mesh = 2m/ (n+1) — 1 G 
[—1, 1] for m = 0, . . . , ri + 1 and n a sufficiently large positive integer. After substituting (f>(ym), rn — 0, . . . , n + 1 
and the approximation to its derivatives in (13 1, we obtain an eigenvalue problem of finite dimension: A^i = cBip, 
for A, B matrices depending only on i?e, a, and n: A is pentadiagonal and B tridiagonal. We solve the eigenvalue 
problem (by means of the inverse power method with adapted shifts) in order to simply get the eigenvalue with the 
largest real part, that is to say, the most unstable one. Precision is improved through extrapolations on the mesh 
size I/in + 1). We have obtained the known results reported by other authors, e.g. Orszag [2], with an analogous 
accuracy. The neutral stability curve, where A^ = 0, is presented in figure [T] In this figure, each point in the Re-a 
plane represents a perturbation of the laminar solution whose stability is decided upon its position: green points are 
stable, red ones unstable and blue ones neutrally stable. We also observe the critical Reynolds number, Recr — 5772.22 
for a = 1.02056, so that if Re < Recr the laminar solution is linearly stable for any value of a. Likewise, for a > 1.1 
the laminar flow is linearly stable for every Re. 



B. Continuation of travelling waves 



Next, we use the above results to look for periodic solutions in time. Due to the translational symmetry of the 
channel in the stream direction (artificial boundaries in ([2|), it is showed in Rand [17| that if we have u{x,y,t) such 
that 

u{x,y,t + T) = u{x,y,t), for all a; e M, ye [-1,1], t^O 
and some T > (that is, u{x,y,t) is T-periodic in time), then it is a rotating (or travelling) wave, i.e. 

u{x,y,t) = u{x - ct,y,0), forc=^. (14) 

Consequently u(x,y,t) is observed as a stationary solution in a system of reference moving at speed c as it was 
introduced in ([5|. The converse is also true, namely every stationary solution of ([5| gives rise to a time-periodic 
solution as is easily verified. This fact allows us to search for periodic solutions in time as functions u{x, y) in a 
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FIG. 2: Bifurcating curve of periodic flows for several discretizations specified as A'^ x M, a — 1.1, and based on Rep and 
Rbq. On each curve based on Rcp there are several corresponding to Hopf bifurcations. They divide the different regions 
of stabihty to superharmonic disturbances, which are also represented in the plot as continuous (stable) and discontinuous 
(unstable) lines. In the Req case, the point labeled Rego represents a real eigenvalue crossing the imaginary axis, meanwhile 
RcQi is a Hopf bifurcation. Likewise, at Rcqq it is attained the minimum Reg. 



Galilean reference at speed 
version of ( 10 ) 



which solve the stationary version of ^ or, in its discretized form, the stationary 



= {I-PkTk)-HUk-PkVk), 



fc = 1,, 



,iV. 



(15) 



Given a fixed a, what we have in (15 1 is a zeros search problem for a system of nonlinear equations of dimension K 
(defined at the end of §111 1. It can be expressed as Hp{Re,c,U) = 0. Solutions of (15) are locally unique for each 
Re, except translations in the stream direction. This is due to the fact that any translation of a rotating wave in the 
stream direction, gives rise to the same wave at a different time instant. Indeed, if u{x,y,Q) is the starting position 
of a rotating wave then, by (14 1, for every 6* g E we have u{x — 0,y,O) = u{x,y,9/c), and thus u{x — 9,y,0) lies in 
the same orbit as u{x, y, 0). In order to achieve uniqueness, we fix one of the coordinates of U, by restricting it to a 
Poincare section 



El = {[/ = {uo, . . .,un) I 3fi(wii) = si} . 



(16) 



(3?(u),S(u) stand for the real and imaginary part of u), for si G M, a fixed value. We mainly set si ~ 0, since this 
choice gives well conditioned systems. If we fix Re in ( 15 ), the number of unknowns, K, is the same as the number of 
equations. We look for its solutions by means of pseudo-Newton's method, in which we factorize the resulting linear 
system using a direct LU decomposition. The first approximation of the Jacobian matrix DHp is implemented by 
finite differences with extrapolation. Every column of this matrix is obtained evaluating Hp in parallel in a Beowulf 
system. Subsequent updates of DHp are carried out by Broyden's 'good' formula. As a consequence, at each pseudo- 
Newton step we only have to apply rank-one updates to the LU factorization of DHp. These improvements mean an 
enormous increase in the speed of computations. 

If we use Re as a continuation parameter, we can trace the one-parameter curve Hp{Re,c,U) — 0. This is 
implemented numerically by pseudo-arclength continuation. To this end, we compute the unit-norm tangent vector 
to the curve. This vector and previously computed points on the curve, are used to predict the next solution point, 
which is finally corrected by pseudo-Newton iterations. 
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FIG. 3: Bifurcating curves of periodic flows for Rep and several values of a specified on each curve. The number of discretization 
points is TV = 22, M = 70. The '•' on each curve represents a change of stability. Joining those points we obtain 4 regions: 
the corresponding solutions are unstable on the right and bottom-left regions and stable on the top region. In the intermediate 
region there are stable and unstable flows, even for a single a, e.g. a = 1.3175. Curves for a = 1.02056 and a = 1.3175 are 
traced in thicker lines. For a = 1.02056 it is attained the critical Reynolds number at Rcp = 5772.22 and for a — 1.3175 it is 
approximately reached a solution at a minimum Rcp = 2939. 



The starting point of those iterations is a numerically integrated periodic solution. We obtain it using the numerical 
integrator (12). The initial condition is taken as a small perturbation of the laminar flow. At the same time we check 
our previous results and those reported in Orszag [2] about Rccr- Namely, taking for example a = 1.02056, Re > Rccr, 
we observe how the laminar flow is not stable, as it evolves to another steady flow which turns out to be T-periodic in 
time for some T. Up to errors of order 0((Ai)^) (those of the time discretization ( 12 1), this T-periodic solution satisfies 



( 15 1 and it is thus a rotating wave. We choose it as initial approximation to a point on the curve Hp{Re, c, U) = 0. 
Given a profile of velocities {u,v) we define its amplitude A, as the distance to the laminar profile (u6,0) in the 



-norm 



A= —\\{U-Ub,v)\\2, 



\\iu,v)\\l = 







[u{x,y)'^ +v{x,yf] dydx. 



(17) 



For a rotating wave as defined in (14), its amplitude does not depend on time since for fixed t: 

p L p\ p L p\ 

I I u{x,y,t)^dydx— / u{x — ct,y,0)^dy dx 
Jo J-i Jo J~i 



L-ct fl 



u{x,y,0) dydi 
1 Jo 



L 1-1 

~ u(i, y, 0)^d?/ di. 

1 
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TABLE I: For some values of a, this table shows the minimum values of Rep and Reg for which exists periodic flow. Calculations 
made for N = 22 and M = 70. The minimum Re attained is marked with Results are in good agreement with those 
reported by Herbert [18] for Rcp. 



min Reg 



min Re„ 



1.1000 
1.2236 
1.3000 
1.3424 
1.3520 
1.3521 
1.3523 
1.3534 
1.5000 



3564.5164 

2845.5884 

2647.6068 

2608.9990 

2607.5519 

2607.5516* 

2607.5520 

2607.5753 

3018.3031 



1.1000 
1.2400 
1.3092 
1.3145 
1.3174 
1.3175 
1.3177 
1.3265 
1.4665 



3797.0331 

3048.0073 

2939.3711 

2939.2069 

2939.0345 

2939.0343* 

2939.0350 

2940.6307 

3526.0725 
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FIG. 4: Real and imaginary part of the two most unstable eigenvalues (Ai and A2) of periodic flows for A'^ — 22, M — 70, 
a = 1.02056 and Rep (cf. figure [3|. Close to the minimum Rep of the amplitude curve, Ai and A2 are unstable and real, and 
give rise to a pair of complex conjugate eigenvalues. The first crossing through the imaginary axis is associated to a single real 
eigenvalue (A2), meanwhile the second one corresponds to the Repi Hopf bifurcation (Ai and A2). Arrows for both graphs, 
point to the direction of increasing amplitude in figure [3] 



In step 1 we apply definition ( 14 ) 
in X we have step 3. 



For step 2 we make the change of variable x 



ct, and because u is L-periodic 



C. Stability of periodic solutions 



We notice that zeros of system ( 15 ) can either correspond to stable or unstable time-periodic solutions. With a 



stable solution it is meant the one for which any small disturbance ultimately decays to zero, whereas if some of those 
disturbances remain permanently away from zero, it is called unstable. 

To decide whether a time-periodic flow u is stable or not, we consider it as a steady solution for its appropriate 
c = L/T and obtain the eigenvalues of its Jacobian matrix. This matrix is computed analytically linearizing (10) 
around u. If every eigenvalue has negative real part, the periodic flow is stable to disturbances of the same wavenumber 
a but, if there is an eigenvalue with positive real part, the solution is unstable. Let us mention that there is always 
a zero eigenvalue which corresponds to the lack of uniqueness of the time-periodic flow due to translations. Setting 
a = 1.1, the bifurcating diagram for the periodic flows in the Re- A plane together with the stability changes are 
represented in figure [2] for both formulations in terms of Rcp and Req. Due to relations (pj) and Q we only need 
to compute travelling waves vP{x,y,t) at speed c for Rcp, since Q gives the associated vP^{x,y,t), which is easily 
checked to be a travelling wave at speed rc for Req, being r = Rep/Req. As well as computing the eigenvalues, we 
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FIG. 5: Real part of the most unstable eigenvalue for periodic flows for a = 1.1 and N x M as specified. For the Rcp case, 
the two crossings of each curve with the imaginary axis corresponds to the first two of figure [2] and are detailed in table [lT| 
as Rcpi and Rep2- In the analogous plot based on Reg, the two crossings of each curve with the imaginary axis corresponds 
to the first two of figure [3] and are also specified in table [ll] as Rego and Regi. Continuous and discontinuous lines refer 
respectively to stable and unstable periodic solutions associated to each point (i?e,3?(A)). Arrows point to the direction of 
increasing amplitude in figure [2] 



confirm the stability of a periodic flow using the numerical integrator (12). A more complete picture of the different 
connections among stable and unstable solutions is given in Casas and Jorba |14j . 

By simply taking a known travelling wave for some a as initial guess and moving slightly a, we can find periodic 
solutions for different values of a. These are shown in figure [3j together with their stability. In addition, in table [T] 
we have computed, for several values of a, the corresponding minimum value of Rcp and Rcq (denoted as Remin{oi)) 
along the amplitude curves (see figure |3|. In turn, Remin{oi), is minimized as a function of a. In this way we obtain 
the absolute minimum Rep and Req for which there exists periodic solution. These minimum values are marked with 
'*' in table|lj For Rep, Herbert HH] obtained the minimum value at Rep = 2934.80 for a = 1.3231 and x M = 4 x 40 
as the spectral spatial discretization for the stream function: this discretization is analogous to the one used in the 
present work. We observe that our value of Rep differ from Herbert's not more than 0.15%. 

Previously, Soibelman and Meiron [T3] found similar bifurcations of travelling waves for a = 1.1, and the critical 
Reynolds number for which there are time-periodic solutions: Rep ~ 2900 for a ~ 1.3, and Reg w 2600. We remark 
that in figure [3] for a = 1.3175 there exists an attracting periodic solution for Rep = 3024 which corresponds to 
Req = 2630. For the case of the laminar flow, we encounter the classical results of Orszag ^2, about the critical 
Reynolds being at Recr = 5772.22 for a — 1.02056. On the one hand, we observe in figure[3]that the bifurcation curve 
of periodic flows reaches the laminar solution at the above mentioned Recr and in addition, the laminar solution is 
checked to be stable when Re < Recr and unstable if Re > Recr- This Hopf bifurcation at Recr is called subcritical, 
because the branch of periodic solutions emanating at it decreases in Re. When the new branch increases in Re, we 
call it supercritical bifurcation. We also notice that for a = 1.02056, it is reached the minimum Reynolds number 
where the transition from stable to unstable laminar flow takes place, what was formerly presented in figure [T] For 
a > 1.1 the curve of periodic solutions does not reach the laminar fiow. This is in agreement with the situation 
shown in figure [l] since for a> 1.1 the laminar flow is linearly stable for every Re. In figure [S] we check as well that 
for a > 0.91 the curve of periodic flows bifurcates subcritically from the laminar fiow, but for a < 0.91 the Hopf 
bifurcation is changed into supercritical. Precisely for a ~ 0.915 it is born a new change of stability on the curve of 
travelling waves at Rep « 6700. The behaviour around this new point is saddle- node bifurcation, analogous to Reqo, 
and will be described in the following subsection. 
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FIG. 6: Speed of the observer c, for parameters and periodic solutions as in figure [2] On both Rcp and -Reg cases, upper 
and lower branches correspond to the respective ones of amplitudes. Unlike Rcp, for the constant flux case the upper branch 
increases with Rbq. Continuous and discontinuous lines refer respectively to stable and unstable periodic solutions associated 
to each c. 



D. Hopf bifurcations 



Now let us discuss the bifurcation diagram shown in figure [3] for a ~ 1.02056. First, the laminar solution becomes 
unstable at the critical value Rccr = 5772.22, due to a Hopf bifurcation that gives rise to a unstable family of periodic 
orbits. This family continues backwards (with respect to Re, i.e. it is subcritical) until Rcp « 4636, where a turning 
point is reached and 3?(A2) (described in figure |4| crosses the imaginary axis. Before arriving at this turning point, 
there is a single eigenvalue, Ai, on the real positive axis, while the remaining ones have negative real part (we ignore 
the eigenvalue at arising from the lack of uniqueness of periodic flows). On traversing through the turning point, 
a real and negative eigenvalue (A2) becomes real positive, so the number of unstable eigenvalues is now two. Shortly 
after that, these two unstable eigenvalues collide and become a conjugate complex pair (still with positive real part), 
and then they cross the imaginary axis for Rcp ~ 4684 producing a new Hopf bifurcation at the point Rcpi on figure [31 
Between Rcpi and Rep2 ~ 6936, the family of periodic orbits is stable to disturbances of the same wavelength. At 
Rep2, there is another Hopf bifurcation produced by a pair of conjugate eigenvalues crossing the imaginary axis. These 
bifurcations persist, as shown in table |llj when Af, N are increased and no new ones seem to appear in this range. 

The case of constant flux is qualitatively different. For Req the bifurcating diagram of periodic solutions has a 
turning point at a minimum value of Reg, which we designate as Rcqq (cf. figure [2]). The lower branch of periodic 
solutions is unstable with only one unstable real eigenvalue and the upper is initially stable, being also real the 
most unstable eigenvalue. On traversing the bifurcating curve towards the upper branch, this real positive eigenvalue 
becomes negative at the turning point. The upper branch is kept stable until a subsequent Hopf bifurcation appears 
at certain value Reqi. For Req > 7000 and a = 1.1, we have detected more Hopf bifurcations which we do not 
consider in this study. However for the range included in figure |2] all periodic flows for Req > Reqi are unstable. 

Pugh and Saffman [TH] pointed out that the null eigenvalue at Reqo has algebraic multiplicity 2 and geometric 
multiplicity 1. We can consider that eigenvalue simple (with algebraic and geometric multiplicity 1) if we ignore 
the constant zero eigenvalue due to a trivial phase shift of the flow in the stream direction. The suppression of this 
trivial null eigenvalue can be made by restricting equations pO| to the closed linear manifold Si defined as a Poincare 
section in (16). According to bifurcation theory, at a simple eigenvalue like this one we have no equilibrium point for 
Req < Reqo and two equilibrium points for Req > Reqo- this situation corresponds to a saddle-node bifurcation 
and no new branches of solutions come out from Reqo. 

In table IH] are shown Rep, Req, the speed of the observer c and the period t of the bifurcated solution corresponding 
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TABLE II; Minimum Reynolds number Rego and the three Hopf bifurcations of periodic flows at Repi, Regi, Rep2, together 
with associated parameters c and r for a = 1.02056, 1.1, M = 70 and several A'^. The values reported in Soibelman and Meiron 
|13| for M = 70 are also included. 
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to the first Hopf bifurcations for several values of N and M — 70. Taking M = 70, Chebyshev modes seem to be 
enough to attain convergence in the results. The values obtained by Soibelman and Meiron il3j are also presented 
for comparison. For a = 1.1 we observe convergence of our results on the different Hopf bifurcations considered as 



N is increased. In all cases there are substantial differences with Soibelman and Meiron s [T^ results, being in more 



agreement for the lowest Re. We remark the slow convergence of the Fourier series to the bifurcation values as N 
is increased. At the same time we have also obtained convergence in the qualitative behaviour: the subcritical or 
supercritical character of all the studied Hopf bifurcations remain unaltered as M, N are increased. 

Formulas ([s]) and Q provide again the correspondence between bifurcation points at Reg and Rcp (cf. table |ll| . 
For instance at Rcpi = 3840.6 for iV = 12, M = 70 and a = 1.1, the periodic solution transformed by these formulas 
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furnish a periodic solution at Req = 3564.5 and c = 0.34490, values in good agreement with Rego reported in table 
[Hj Likewise, the transformed periodic solution for N = 25, a = 1.1 at Reqi = 5249.6 gives rise to Rep = 7565.5 and 
c = 0.28192, again in good precision with respect to Rep2- 

The different stability changes marked as a blue dot in figure [3j roughly divide the Rcp-A plane in four regions. On 
the right and bottom- left regions the corresponding periodic solutions are unstable (filled with red points), meanwhile 
they are stable on the top region (only green points). In the intermediate region there are both stable and unstable 
flows, even for a single a, e.g. a — 1.3175. Through this classification, given a periodic flow with its associated 
(Rep, A), we can deduce its stability, independently on a in some cases. On traversing the Rcpi blue curve in the 
direction of increasing amplitudes, up to the relative maximum on that curve attained at Rcpi ~ 9909 for a ~ 0.814, 
we encounter saddle-node bifurcations for a G [0.88, 0.915] at a relative maximum on each bifurcation curve. For 
a G [0.714, 0.88], the former relative maximum disappears and the saddle-node bifurcation turns into a Hopf one. The 
rest of the Rcpi curve is made up of Hopf bifurcations for the studied values a ^ 1.7. The minimum Repi « 3024 is 
reached precisely for a ~ 1.3175, where the minimum periodic flow was found in § IV C The Rep2 blue curve is only 



constituted of Hopf bifurcations for the a e [0.74, 1.3175] considered. The minimum Rep2 ~ 6936 is achieved again 
for the critical a « 1.02056. 

The maximum growth rate (real part of the most unstable eigenvalue) for each periodic flow is presented in figure [5] 
for the same parameters as in figure [2j For the most unstable eigenvalue A, 3?(A) crosses the imaginary axis twice, on 
the values Repi and Rep2 for Rep and on Rego and Rbqi for Req. Those diagrams represent the degree of instability 
of each flow. Comparing to figure [2| we observe that at the same Re on the upper branch of amplitudes, periodic 
solutions based on Rcq are more unstable than the associated ones based on Rep. On the other hand, on the lower 
branch of amplitudes, at the same Re, both curves of 3fi(A) visually coincides for Re > 5000. This behaviour is also 
refiected in figures [2] and [6j In this last figure we present qualitatively different curves for the speed c in Rep and Reg 
cases. In the first case, c is decreasing in both branches of solutions in figure [2] However, for Req the shape of the 
c-curve is similar as the A-curve in figure [2] 



V. QUASI-PERIODIC SOLUTIONS 



In this section, we study the quasi-periodic flows that appear at the Hopf bifurcations of rotating waves shown in 
§ |IV[ They are found as time-periodic orbits in an appropriate Galilean reference, which simplifies enormously their 
search. Those time-periodic orbits are obtained as fixed points of a Poincare section, by means of a pseudo-Newton 
method. We have traversed bifurcating branches of quasi-periodic solutions for the Hopf bifurcations at Repi, Rep2 
and Reqi defined in §IV We have obtained different qualitative results than the ones reported in Soibelman and 



1.1 and Rep they found that the bifurcation at Repi to quasi-periodic solutions is subcritical: 



Meiron [13]. For a 

they obtained quasi-periodic solutions for Rep before the bifurcation point. In consequence, close to that point, those 
bifurcated flows are stable. In the present study, by increasing the number of Fourier modes N, we have achieved a 
supercritical Hopf bifurcation at Repi: the bifurcating quasi-periodic flows are located for Rep after the bifurcation 
point and therefore close to it they are unstable. This is treated in § VB[ For the second Hopf bifurcation at Rep2, in 
agreement with Soibelman and Meiron |13j . the quasi-periodic orbits are found for Re 
More details are given in 



greater than the bifurcation 
VC The behaviour at Reqi (considered in § VD) is analogous to that of Rep2- 
> Reqi large enough we have detected another Hopf bifurcation to tori with 3 basic frequencies. 



point 

However for Req 

The stability of quasi-periodic flows to superharmonic disturbances is estimated by means of the linear part of the 
Poincare map and also with a full numerical simulation of the fluid. 



A. Reduction to periodic and numerical procedures 



We use again the spatio-temporal symmetry of our system, due to the artificial boundaries of the channel. Consid- 
ering this symmetry in Rand [17] it is proved that every solution u{x, y, t) that lies on an isolated invariant 2-torus 
(a quasi-periodic solution), not asymptotic to a rotating wave, is a modulated wave, that is to say, there exists t > 
and € M such that 

u{x, y, nr + t) = u{x — ncf), y, t) for every n G 1j. (18) 

Hence, this kind of wave has the property that, may be viewed as a r-periodic wave in time, in a frame of reference 
moving at speed c = {pL + (j))/T, for any integer p. In effect, defining x = x — at for that value of c, and u{x, y, t) = 
u{x + at, y, t) as the velocity in the moving frame of reference at speed c, it turns out that 

u{x,y,T) = u{x + cT,y,T) = u{x + pL + (l),y,T) = u{x,y,Q) = u{x,y,0). (19) 
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FIG. 7: Two representations of the solution vector U{t) of a quasi-periodic flow projected on the plane of two selected coordinates 
namely, Af - 1 + (2iV - 3)(A/ - 2) + 3, and M - 1 + 2(M - 2) + 5, for Rep = 3865, a = 1.1, iV = 8, Af = 70. The range of 
values for (a) is [-0.002,0.002] x [-0.001,0.001] and [0.00037,0.00129] x [0.000338,0.000385] for (b). In (a) a dot is plotted 
each At — 0.02 time units, meanwhile in (b) only for t such that U{t) £ Ei. The red closed curve in (b) is obtained integrating 
repeatedly for t G [0, r]. This curve is also represented as a reference on the right centre of (a). 



In steps 1, 4 we use the previous definition of u. Substituting the previously defined c we obtain step 2, and because u 
is L-periodic in x and a modulated wave, one gets step 3. Consequently we have proved that u{S:, y, t) is a r-periodic 
function of t. 



In order to look for periodic flows satisfying (191 in a Galilean reference at speed c, we make use of the Poincare 
section Ei defined in (16). In this case, we only consider points on Si when they cross the section in a particular 
direction as time increases, namely, from si < to si > 0. Likewise we define the associated Poincare map Pc '■ 
Si — > El as follows: starting from an initial condition U = U{0) G Ei we integrate (10 1 for fixed parameters Re, a 
and c, until a time such that U{tc) & Ei {U{tc) represents the evolution of U{0) in a Galilean reference at speed 
c) for the nc-th time (i.e. after ric crosses with Ei), where ric is a positive integer which represents the minimum 
number of times needed for the flow to return close to the initial point U{0) (the meaning of 'close' will be specified 
in §VB). We then set Pc(J7(0)) = U{tc)- In this way we have reduced the search of quasi-periodic flows to a zeros 



flnding problem for the map Hg defined as 



= Hg{Re, c, U) = Pc{U) -U = U{Q - U(0). 



(20) 



From (pi we have that, if Pc{U{0)) = U{t) = {7(0), for some r, then as well PciV{0)) = V{t) = V{0) for 
V{0) = P^U{0)) and k any positive integer. Indeed, we can express ( [l9| as u{x,y,T + ti) = u{x,y,ti) for every 
ti. From here, since 1^(0) = U{ti) for some ti, we immediately obtain P^V{0)) = V{t) — V{0). Therefore we can 
generate different points on the same orbit as a solution of (20 1. We avoid this lack of uniqueness by restricting 
Pc : TiiD E2 — y El, for E2 a Poincare section (analogous to Ei) defined by 



T,2 = {U = (wo, ■ • ■ , un) I S" = 0} , 



(21) 



where we set S = 3?(uAr.M/2-i) ~ S2, for S2 G K a suitable quantity. For Re close to the studied Hopf bifurcations, we 
have chosen S2 = ^{^^ ^1/2-1"^^ with £ Ei the travelling wave at the exact Re where the bifurcation takes place. 
The reason for this choice is merely to preserve continuity in the amplitude diagrams described next. 

In order to trace the curve Hg{Re, c, U) = 0, we utilize the same continuation method as for system Hp{Re, c,U) ~ 
in (15), differing essentially in the definition of the equation to vanish: for periodic flows the computations are much 
simpler and faster than for quasi-periodic ones. The solution of (20) needs an initial guess, which is obtained as 
described in the following subsections. Once we have a quasi-periodic flow such that Hq{Re, c,U) = 0, we measure 
its amplitude A, as in the case of periodic flows (cf. ([l7|). If u{x,y,t) is a modulated wave, using (18) and the 
L-periodicity we have 



u{x, y, r)^ dy dx = 



1 i<L 

u{x - (j),y,O fdydx = 
1 Jo 



u{x, y, 0)^ dy da:. 



(22) 
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FIG. 8: Bifurcated branches of quasi-periodic flows at the supercritical Hopf bifurcation of periodic flows at Rcpi. Each curve 
represents A as a function of Rep and has specified N x M. Both branches of periodic (in lighter colors) and quasi-periodic 
flows are presented. Calculations are shown for a = 1.1 and At = 0.02. The correspond to the Hopf bifurcation at Rcpi 
presented in figure [3] For the range shown, the bifurcating branch consists of unstable quasi-periodic orbits (dotted lines) , 
whose amplitude decreases with Rep. 



Since we have numerically checked that A is not constant for modulated waves, we conclude from (22) that it is a 
T-periodic function of t. This is not so for the rotating waves of § IV for which A is constant in time. In the case 



of a quasi-periodic flow U{t), with the purpose of considering a concrete value for the amplitude, we evaluate A(t), 
at t such that U{t) e Ei H T.2. This is simply a representative and easy to compute value for A(t), as we cannot 
obtain a single value for the amplitude of this class of flows. We use it to trace the continuation curve: it provides 
the distance to the laminar solution at some time instant. In the same way, c can be considered as a representative 
and time independent value for every quasi-periodic flow, so that we can as well use it to trace continuation curves. 



B. Hopf bifurcation at Repi 



For Rcp < Rcpi and a « 1 the corresponding time-periodic flow in figure^ is unstable, but its temporal evolution 
ultimately decays to the laminar flow. The results in Soibelman and Meiron [13] point out the existence of a subcritical 
Hopf bifurcation at Repi: they use the vorticity equation and only consider N ^ 2 Fourier modes. According to 
bifurcation theory (see Marsden and McCracken [20]), this means that the bifurcating quasi-periodic flows are locally 
stable. Following this result we tried to find quasi-periodic flows in the subcritical region, but with no success: we 
were not able to detect a quasi-periodic attracting solution for Rcp < Rcpi and TV ^ 3. In consequence we direct the 
search of quasi-periodic flows to the supercritical region, i.e. for Rep > Rcpi. 

It is also known from bifurcation theory that, in the case of a supercritical Hopf bifurcation in which flxed points 
before it are unstable and after it stable, the branch of periodic solutions that emanates at the bifurcation point is 
locally unstable. That should be the situation for Rcp > Rcpi and thus the bifurcating quasi-periodic flows are locally 
unstable and therefore hard to locate by direct numerical integration, because the evolution of the fluid close to them 
does not remains near as time evolves. 



On the other hand, pseudo-Newton's method applied to solve (pO|) does not distinguish between stable or unstable 

it is carried out as follows. First we 
Close to the 



flows. However, the difficult task is the search of a good starting guess for ( 20 ) 

consider an unstable time-periodic solution Uf(t), i.e. a fixed point of (10 1, for certain Rei < Rcpi 



Hopf bifurcation Rcpi in the unstable region Rcp < Rcpi, the linearization of (10) around Uf has just a couple of 



complex conjugate eigenvalues, A = ± iX^ with A^ > (cf. § IV). If w = ± iwi is the associated eigenvector, we 
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FIG. 9: Different curves related to figure[8]for N x M points as specified, (a) r represents the period in time for a quasi-periodic 
flow when the observer sees it as periodic, (b) The appropriate value of c that converts a quasi-periodic flow in periodic, together 
with the corresponding curve for periodic flows in lighter colors (cf. figure |6|. The '*' correspond to the Hopf bifurcations at 
Rspi presented in figure [2] 




FIG. 10: Two representations of the solution vector U{t) of an (almost resonant) quasi-periodic flow projected on the plane of the 
same coordinates as figure[7| for Re^ = 7000, a = 1.02056, iV = 22, M = 70. The range of values for (a) is [-0.00185, 0.00185] x 
[-0.004,0.004] and [0.00147,0.00164] x [-0.002495,-0.002385] for (b). In (a) a dot is plotted each At = 0.01 time units. In 
(b) a dot is plotted only for t such that U{t) £ Ei. The curve in (b) is also represented as a reference on the lower right corner 
of (a) as a small green line. 



choose V e {wr,Wi), i.e. v is in the plane of the most unstable directions, along which the flow escapes in the fastest 
fashion from C/f . Now we change to Re2 > Rcpi close to the bifurcation and select C2 near to the one associated 
with the corresponding periodic orbit U2 at Re2, i.e. Hp(Re2,C2, J/f) = 0. For Re2,C2 and jrj <C 1 a small constant, 
we integrate numerically the initial condition U = Uf -\- rv, until a time when it become as closest as possible to a 
solution of Hq{Re2, C2,U) = 0. This first approximation of C2 is optimized by using minimization algorithms. In this 
case we also observe a value t of the return time of Pc, neighbouring to 27r/|Im(A)|, taking A as the unique purely 
imaginary eigenvalue (together with its conjugate) at the Hopf bifurcation Rcpi. 

The initial guess obtained in this way is improved through pseudo-Newton's iterations applied to the function Hq, 
to finally obtain a first unstable modulated wave for Rcp > Rcpi . In figure [7| we present numerical evidence that 
solutions U{t) of (20) lie in a 2-torus and are unstable in time for supercritical Rcp. In a) we observe how the 
trajectory, projected on a plane of two arbitrary coordinates of U{t), fills densely a 2-torus. On the outer red curve 
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FIG. 11: Bifurcated branches of quasi-periodic flows at the supercritical Hopf bifurcation of periodic flows at Rep2- Each curve 
represents y4 as a function of Rep and has specified A'' x M. Both branches of periodic (in lighter colors) and quasi-periodic 
flows are presented (the green color almost completely conceals the red one). Calculations are shown for a = 1.02056 and 
At — 0.02. The correspond to the Rcpi, Rep2 Hopf bifurcations presented in figure [s] For the range shown, the bifurcating 
branch consist of quasi-periodic orbits, stable to disturbances of the same wavenumber. 



of b), it is simply plotted U{t) when U(t) G Si with the following adaptation: ii t > t, being r the period of U{t) as 
a modulated wave, we plot U{s) for s such that 0^s = t — pT<T and p ~ [t/r], {[■] stands for the integer part) i.e. 
we treat U{t) as if it were exactly r-periodic in order to avoid its unstability. We note that, as we are on the Poincare 
section Si, the 2-torus in a) is reduced to a closed curve in b), which corresponds to the unstable quasi-periodic flow, 
seen as if it were unperturbed by numerical errors. Likewise in b) we let the flow evolve for long time and plot again 
U{t) (in green) when U{t) € Si but suppressing the previous time adaptation. Here we can check that the flow is 
unstable, because it moves away from the outer closed curve and falls to the time-periodic and stable regime: a simple 
point in the centre of the figure. 

Once we have obtained the first solution of (20) by the pseudo- Newton's method, we use continuation methods to 
traverse the bifurcating branch of quasi-periodic flows parametrized by Rcp. In figure l8| we plot the amplitude A for 
each quasi-periodic solution as a function of Rep. It seems that we have achieved both qualitative and quantitative 
convergence because we obtain a similar graph in increasing the values of N, M. It looks also clear from this plot 
that the Hopf bifurcation is supercritical so, at least locally, solutions on the bifurcating branch are unstable. The 
analysis of the stability of a quasi-periodic solution is done by means of the eigenvalues of the linear part of Pc at 
a fixed point U{t) such that Pc{U{Q)) — U{0). This computation, analogously to the matrix DHq, is obtained by 
extrapolated finite differences. In the range of Rep presented in figure [Sj the quasi-periodic solutions are unstable. 



Furthermore, unlike the bifurcation at Rep2, solutions at Repi present certain symmetry: Uk{y,t) (defined in §111) is 
an even or odd function of y according to the parity of k. 

We also observe in figure|9|^a) for different values of N, an indicator of the numerical effort involved in the evaluation 
of the map Pc- The big slope of the Rep-r curve shows the high computational cost involved as Rep is slightly increased. 
For Rcp e [3840,3872], a = 1.1 and At = 0.02 the time needed to return to Si is r S [7000, 10700]. Likewise it has 
been necessary to adapt dynamically the minimum number of times, ric (defined just before system (20l), that the 
solution crosses Si before it returns to the starting point. For Rep close to Rcpi, Uc starts at 1 and it is incremented 
by 1 when Rep varies approximately in just 2-3 units. The way we modify is described next. If we change slightly 
Rep the return time r should be varied in accordance. Thus, we impose that the new value obtained for r satisfy 
|r — Toj < £, for some tolerance e and Tq the former return time. If this condition is not fulfilled we increase or decrease 
ric by one unit, until it is satisfied, or we decrease Rep if necessary. 

In figure l9|b) is represented c for the same range of Rcp. We observe for c, nearby values as their counterpart 
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FIG. 12: Different curves related to figure [Tl] using Rep in the abscissa axis for N x M points as specified, (a) and (b) are as 
in figure|9] In (b), in lighter colors, it is also drawn the periodic flows around Rep2 (marked with '*'). In this case, the green 
color almost completely conceals the red one. 
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FIG. 13: Analogous to flgure [TT| based on Regi, using the specifled a and N x AI — 22 x 70, At = 0.01. The bifurcated 
branches of quasi-periodic solutions have a change of stability at another Hopf bifurcation. Continuous and discontinuous lines 
represent stable and unstable flows respectively. 



periodic flows (cf. figure |6|, but in this case the curve has a large decreasing slope. We emphasize that the c graph 
shows a bifurcation diagram of periodic and quasi-periodic flows, which is time independent in contrast to the Re- A 
plot in figure [8] 




C. Bifurcation at Rep2 



As we presented in figure |3j for Re > Rep2 the corresponding periodic flow is unstable, so tfie evolution of ( 10 ) 
from such a flow as initial condition, drives the fluid away from it. By following the temporal evolution of this flow, 
we observe that the fluid seems to fall in a regular regime, which flnally proves to be a quasi-periodic attracting 
solution. This is checked in figure 10 where we plot the projection of the solution vector U (t) over the plane of the 
same two coordinates as in figure [tj Each point in figure 10 'a) corresponds to the value of the specified coordinates 
at a time instant. As we can observe, the trajectory appears to fill densely the projection of a 2-torus. Plotting the 
same coordinates as above, but only when U{t) e Ei, we see in figure ^h) a closed curve, which seems again to 
confirm that the flow lives in a 2-torus. 

Let us denote as U^{t) = (ug, . . . , u^){t), a time instant of the flow in the attracting 2-torus. We need to approximate 
the value of c° which better makes appear as a periodic flow. That c° exists according to (18|-(19). It can be 
estimated as (cf. Rand [T7] ) 



27r n{t) 
c = — lim , 

k i->oo t 



for n{t) = 



arg(Conj«„(i))) 



2n 



where k is the number of peaks of the wave C/" for x G [0, L] and we have used the midpoint of the channel for 
m = M/2. With t large enough, 2Trn{t)/kt gives an approximation of c" which we optimize using minimization. Next 
we use U'^, c° as the initial condition for pseudo-Newton's method applied to ( 20 ) to conflrm that our attracting 2-torus 
C/° is in fact a modulated wave. For a flxed a, once we have a first point {Rcp, c", JJ") which satisfies (20), taking Rcp 



as a continuation parameter, we can trace the curve of quasi-periodic fiows in the Re~A plane by pseudo-arclength 



numerical continuation, applied to Hq as in § V B 



As we observe in figure |11[ there appears a branch of quasi-periodic solutions which bifurcates supercritically 
from the curve of periodic flows. We may again have zeros of Hq(Re,c,U) that can either correspond to stable 
or unstable quasi-periodic solutions. Since on crossing the bifurcation point Rep2, the periodic orbits change from 
stable to unstable, the branch of bifurcating quasi-periodic flows are locally stable to two-dimensional superharmonic 
disturbances. By means of the eigenvalues of the Jacobian matrix dPc/du we also compute the stability of the obtained 
quasi-periodic solutions. For a — 1.02056 and the range of Rcp e [7000, 13000] studied, all quasi-periodic flows found 
are stable to perturbations of the same wavelength and the situation is kept when N, M are increased. In flgure [12] 
we present the curves of frequencies c, t which deflne the different modulated waves. Again the Re-c graph shows a 
time independent bifurcation diagram. 
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FIG. 15: (a) An stable 2-torus on Si for Reg = 8635, a = 1.02056, iV x = 22 x 70, At = 0.009. The integration time 
on the figure is about 1,222,000 time units, (b) Analog of (a) for Req — 8640. We observe in this case that the initial 
2-torus is unstable and is attracted by a 3-torus. The integration time on the figure is about 419,000 time units, (c) The 
same as (b), but the initial condition is an unstable 2-torus for Reg — 9750 which is also attracted by a 3-torus. The 
integration time is about 1,005,000 time units, (d) The initial condition is an unstable 2-torus for Rcq — 10500 and a — 1.10 
which is again attracted by a 3-torus. The integration time is about 1,002,000 time units. The respective ranges in the four 
figures are [0.0019, 0.00445] x [-0.00965, -0.0053], [0.0018, 0.0045] x [-0.0098, -0.0052], [0.0013, 0.0058] x [-0.0126, -0.005] and 
[0.0016,0.00535] x [-0.0147,-0.0078]. 



D. Bifurcation at Regi 



In the case of constant flux the bifurcation diagram of periodic flows is qualitatively different to that of constant 
pressure, as can be verified in figure[2] For Reg and the values of a considered (0.9, 1.02056, 1.1 and 1.15) there is a 
change of stability at the minimum Reynolds Rego of the amplitude curves, but no new bifurcations are born there. 
The first Hopf bifurcation occurs at the point labeled Rcqi in figure[2j We can summarize that the qualitative picture 
of amplitudes of quasi-periodic solutions emanating from Rbqi is analogous to that of Rep2. The main differences are 
basically quantitative, because Reqi < Rep2 (cf. figures 
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and 13 1 



For a = 0.9, 1.02056, 1.1 and 1.15, we compute the quasi-periodic flows that bifurcate from Rcqi, following the 
same steps of § V C In figure 13 we plot their amplitudes, and again, as in the case of Rep2 we observe a supercritical 
bifurcation. The associated frequencies c, r, to the modulated waves are presented in figure 14 We remark the 
analogies between bifurcation diagrams of A and c in respective figures 13 and [l4]^b), the latter one being time 
independent. The quasi-periodic solutions found from Reqi and a — 1.02056, are stable for Reqi < Reg < 8640. At 
Reg 8640 the branch of quasi-periodic solutions loses stability at a new Hopf bifurcation, giving rise to a family 
of attracting tori of 3 frequencies. Numerical evidence of this bifurcation is shown in figure |15[a) where the same 
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FIG. 16: Stream lines (left) and levels of vorticity (right) for a unstable quasi-periodic flow at Req — 10500, a — 1.1, A'^ = 
22, M = 70, At = 0.009. Surface levels are plotted for ^ x ^ L, —1 ^ y ^ 1. The depicted time instants from top to down 
are 0, r/6, 2r/6, 3r/6, 4r/6, and 5r/6 for the associated orbit period r — 10.24. Ranges for levels of stream lines and vorticity 
are [—0.55, 0.80] and [—3, 3] respectively. They are represented by a common scale of colors to the right of each column. 



coordinates of U{t) as in figure [t] are plotted on Ei for Req = 8635, yielding an apparently perfect closed curve, after 
a considerably long time of integration: it is a stable quasi-periodic solution. On the contrary, in a similar plot for 
Req — 8640, figure 15 h) shows an unstable quasi-periodic flow, which is attracted by a 3-torus. The new frequency 
of this attracting solution is verified using in turn the Poincare section E2. We have carried this out by plotting the 
same two selected coordinates of U{t) only when the flow crosses Ei if, in addition, it is approximately on E2. We 
obtain in this way what seems to be a closed curve (see Casas and Jorba [T4l for a plot) as may be expected for a 
3-torus. We can observe two other more involved 3-torus in figure 15 'c) and (d). 
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Finally in figure [THj stream lines and levels of vorticity of an unstable quasi-periodic flow for Req = 10500, a = 1.1, 
X M = 22 X 70 and At = 0.009, are presented at six equidistant values of time in the interval [0, r] for r = 10.24. 

We observe (as expected) larger vorticity close to the walls than in the channel centre, which in addition can be 

confirmed on the stream lines figures. 



VI. CONCLUSIONS 



In this work we have studied some bifurcations of plane Poiseuille flow. We have reproduced results of other authors 
and obtained similar qualitative results about the Hopf bifurcations, in what concerns to their number and location. 



The main quantitative differences between Soibelman and Meiron s [13] computations and ours are due to the larger 



resolution we have used, together with the distinct formulations implemented of the Navier-Stokes equations. The 
important qualitative difference is the kind of bifurcation found at Rcpi: in their computations this bifurcation is 
subcritical, but improving the precision of the numerical approach we obtain that it is supercritical. Then, the 
bifurcating quasi-periodic orbits are unstable. This has also been conflrmed by numerical simulations. 

In the case of the bifurcation at Rcpi, because the lengthy time integrations, we have only been able to move away a 
few tens from Rcpi . The further we advanced in i?ep, the greater are the numerical difficulties we encounter to track the 
bifurcating branch of quasi-periodic flows, due to long time integrations. It is also worth to mention the complications 
derived from their instability. Close to Rcpi and with the discretization employed {N = 8, M — 70, At — 0.02), it 
seems that we have achieved both qualitative and quantitative convergence. By observing figure [8] we can conjecture 
that, for the range of a G [1, 1.1] considered, the minimum Re « 2900 attained with travelling waves is not lowered 
by quasi-periodic fiows. This question still remains open for two-dimensional flows, although Ehrenstein and Koch 
[S] solved the gap between experiments and numerical results in the case of three-dimensional flows. However, in the 
present work for a = 0.89, Rcpi ~ 7250 we have localized a subcritical branch of stable quasi-periodic orbits. They 
are difficult to follow by means of the approach described in § |VB[ because the time needed to evaluate the Poincare 
map is r « 15, 000 time units. It remains open whether that family could reduce the minimum Re 2900 of periodic 
flows to Re « 1000, where transition has been observed experimentally. 

For Rcp > Rep2 the quasi-periodic flows encountered are attracting and the integration time is of the order of 
tens, so in this case the computational cost is drastically reduced compared with the bifurcation at Repi. The 
range of Rep obtained for attracting quasi-periodic flows moves now to several thousands. However, in spite of 
keeping qualitative convergence, the use of larger Reynolds numbers makes necessary an increase in precision to get, 
furthermore, quantitative convergence. An analogous qualitative picture is found at the quasi-periodic flows which 
bifurcate from Reqi. In this case the quasi-periodic solutions quickly loses stability and we have also obtained another 
Hopf bifurcation to a family of tori with 3 basic frequencies. We could say that dynamics are richer for Req than 
Rep (see Casas and Jorba il4j), because bifurcations and different vortical states appear for lower Req than the 
counterpart Rep. 

As future work, it would be of interest to analyse the stability to disturbances with different wavenumber a or even 
to 3-dimensional perturbations. The connections of the different families of solutions is also of great relevance, or 
even the discovering of new vortical states which could approach more the transition to turbulence. Likewise, due to 
nonnormality in the Navier-Stokes system, the sensitivity of eigenvalues to perturbations is an important issue that 
can be analyzed by means of the pseudospectra. This study would give a measure of the reliability of the spectrum 
obtained in linearizing ( 10 1 around periodic flows, mainly for high values of Re. In this work, the stability according 
to eigenvalues is coherent with our direct numerical simulations ( 12 ) of the different flows. 
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